diary cpresult
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Figure A.1: peaceful power transition, quick health deteriorate
clear

T=1000;
b=10;
del=0.7;
R=10;
r = 0.1;
l = 1;
L = l*0.001;
h = 10;
H = h*0.001;
ph=0.5;
w=0.05;
eta = 0.7;
S_bar = 0.01;

p = [0:0.01:1 ones(1,900)];
s=0:0.001:1;
q=1-p;


k=zeros(1,T);
for i=1:T
    AA=0;
    for j=i+1:T               
        A=del^(j-i)*prod(q(i+1:j));
        AA=AA+A;
    end
    k(i) = r+r*AA;
end


Vm_4=zeros(length(s),T);
Vm_3=zeros(length(s),T);
Vm_2=zeros(length(s),T);
Vm_1=zeros(length(s),T);
Mcut = zeros(T,1);
Vc_4=zeros(length(s),T);
Vc_3=zeros(length(s),T);
Vc_2=zeros(length(s),T);
Vc_1=zeros(length(s),T);
Ccut = zeros(T,1);
Eff=zeros(length(s),T);
for t=1000:-1:1
 % ruler's decision at stage 4
  % if keeps,
 keep=zeros(1,length(s));
 for jj=1:length(s)     
     s4 = s(jj);
     if t<1000
         keep(jj) = r+del*(p(t+1)*min(s4+w,1)*eta*R+(1-p(t+1))*Vm_2(jj,t+1));
     else
         keep(jj) = r+del*(p(t+1)*min(s4+w,1)*eta*R);
     end
 end
  % if strips
 strip = k(t);  
 dif=keep-strip;
 % check uniquness
 nc=0;
 po=zeros(1,length(s));
 for vv=2:length(dif)
     if dif(vv-1)>0 && dif(vv)<=0
         nc=nc+1;
         po(nc)=vv;
     end
 end
 if nc==1
     Mcut(t) = s(po(nc));
 elseif nc==0
     Mcut(t) = 1.1;
 end
 for vv=1:length(s)
     if s(vv)<Mcut(t)
         Vm_4(vv,t)=keep(vv);
     else
         Vm_4(vv,t)=strip;
     end
 end
% cp's payoff at stage 4
 for jj=1:length(s)
    s4=s(jj);
     if s4>Mcut(t)
         Vc_4(jj,t)=-b;
     else
         if t<1000
             Vc_4(jj,t)=del*(p(t+1)*(min(s4+w,1)*(R+b)-b)+(1-p(t+1))*Vc_2(jj,t+1));
         else
             Vc_4(jj,t)=del*(p(t+1)*(min(s4+w,1)*(R+b)-b));
         end
     end
 end
 
%cp's decision at stage 3
lo = zeros(1,length(s));
cl = zeros(1,length(s));
for jj=1:length(s)
 % challenge
 s3=s(jj);
 cl(jj) = min(s3,1)*(R+b)-b;
 % loyal
 if t<1000
     lo(jj) = Vc_4(jj,t);
 else
     lo(jj) = del*(p(t+1)*(min(s3+w,1)*(R+b)-b));
 end
end
 cdif=lo-cl;
  % check uniquness
 nc=0;
 po=zeros(1,length(s));
 for vv=2:length(cdif)
     if cdif(vv-1)>0 && cdif(vv)<=0
         nc=nc+1;
         po(nc)=vv;
     end
 end
 if nc==1
     Ccut(t) = s(po(nc));
 elseif nc==0
     Ccut(t) = 1.1;
 else
     disp('multi-cut for cp');
     t
     test=po(po>0);
     length(test)
     for qq=1:length(test)
         [lo(test(qq)) cl(test(qq))];
     end
     Ccut(t) = s(po(1));
 end
 for vv=1:length(s)
     if s(vv)<=Ccut(t)
         Vc_3(vv,t)=lo(vv);
     else
         Vc_3(vv,t)=cl(vv);
     end
 end
 % m's payoff at stage 3
 for vv=1:length(s)
     if s(vv)<=Ccut(t)
         Vm_3(vv,t)=Vm_4(vv,t);
     else
         Vm_3(vv,t)=-(b+k(t))*min(s(vv),1)+k(t);
     end 
 end
 
 % cp's decision at stage 2. 
  for jj=1:length(s)
      s2=s(jj);
     % e=h
      s3_h=s2+H;
      p3_h=min(jj+h,1000);
      s3_l=s2+L;
      p3_l=min(jj+l,1000);
      Eh = ph*Vc_3(p3_h,t)+(1-ph)*Vc_3(p3_l,t);
     % e=l
      El = Vc_3(p3_l,t);
     % decision 
      if Eh>=El
          Vc_2(jj,t)=Eh;
          Eff(jj,t)=1;
      else
          Vc_2(jj,t)=El;
          Eff(jj,t)=0;
      end
  end
  
  % ruler's payoff at stage 2
   for vv=1:length(s)
       
     if Eff(vv,t)==1
         Vm_2(vv,t)=ph*Vm_3(min(vv+h,1000),t)+(1-ph)*Vm_3(min(vv+l,1000),t);
     else
         Vm_2(vv,t)=Vm_3(min(vv+l,1000),t);
     end 
   end
end  
% Find critial points 
  % t_hat
AA = Mcut-Ccut;
BB = zeros(length(AA),1);
AA(AA<=0)=-1;
AA(AA>0)=1;
t_hat = find(AA>0,1)

  % t_m
t_m = find(Mcut>1,1) 

% graph
 % conflict thresholds
CC = Mcut;
d=(Mcut(t_hat)-Mcut(1))/(t_hat-1);
V1=Mcut(1):d:Mcut(t_hat);
CC(1:t_hat)=V1;
for i=1:length(Mcut) 
    if Mcut(i)-Ccut(i)>0.045 || Mcut(i)>=1
        CC(i)=Ccut(i);
    end
end

  % power increase
PI = (ph*(H-L)-L)*(1:T)'+S_bar;
QQ = CC-PI;
End = find(QQ<=0,1);

plot(1:End-1,CC(1:End-1));
hold on
plot(1:t_hat,CC(1:t_hat),'r')
plot(1:End-1,PI(1:End-1),'g');
hold off

%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Figure A.2: moderate health deterioate, peacefull power transition.
 
clear

T=1000;
b=10;
del=0.7;
R=10;
r = 0.1;
l = 1;
L = l*0.001;
h = 10;
H = h*0.001;
ph=0.5;
w=0.05;
eta = 0.7;
S_bar=0.29;


p = [0:0.005:1 ones(1,800)];
s=0:0.001:1;
q=1-p;


k=zeros(1,T);
for i=1:T
    AA=0;
    for j=i+1:T               
        A=del^(j-i)*prod(q(i+1:j));
        AA=AA+A;
    end
    k(i) = r+r*AA;
end


Vm_4=zeros(length(s),T);
Vm_3=zeros(length(s),T);
Vm_2=zeros(length(s),T);
Vm_1=zeros(length(s),T);
Mcut = zeros(T,1);
Vc_4=zeros(length(s),T);
Vc_3=zeros(length(s),T);
Vc_2=zeros(length(s),T);
Vc_1=zeros(length(s),T);
Ccut = zeros(T,1);
Eff=zeros(length(s),T);
for t=1000:-1:1
 % ruler's decision at stage 4
  % if keeps,
 keep=zeros(1,length(s));
 for jj=1:length(s)     
     s4 = s(jj);
     if t<1000
         keep(jj) = r+del*(p(t+1)*min(s4+w,1)*eta*R+(1-p(t+1))*Vm_2(jj,t+1));
     else
         keep(jj) = r+del*(p(t+1)*min(s4+w,1)*eta*R);
     end
 end
 strip = k(t);  
 dif=keep-strip;
 % check uniquness
 nc=0;
 po=zeros(1,length(s));
 for vv=2:length(dif)
     if dif(vv-1)>0 && dif(vv)<=0
         nc=nc+1;
         po(nc)=vv;
     end
 end
 if nc==1
     Mcut(t) = s(po(nc));
 elseif nc==0
     Mcut(t) = 1.1;
 else
     disp('multi-cut for M');
     t
 end
 for vv=1:length(s)
     if s(vv)<Mcut(t)
         Vm_4(vv,t)=keep(vv);
     else
         Vm_4(vv,t)=strip;
     end
 end
% cp's payoff at stage 4
for jj=1:length(s)
    s4=s(jj);
     if s4>Mcut(t)
         Vc_4(jj,t)=-b;
     else
         if t<1000
             Vc_4(jj,t)=del*(p(t+1)*(min(s4+w,1)*(R+b)-b)+(1-p(t+1))*Vc_2(jj,t+1));
         else
             Vc_4(jj,t)=del*(p(t+1)*(min(s4+w,1)*(R+b)-b));
         end
     end
 end

%cp's decision at stage 3
lo = zeros(1,length(s));
cl = zeros(1,length(s));
for jj=1:length(s)
 % challenge
 s3=s(jj);
 cl(jj) = min(s3,1)*(R+b)-b;
 % loyal
 if t<1000
     lo(jj) = Vc_4(jj,t);
 else
     lo(jj) = del*(p(t+1)*(min(s3+w,1)*(R+b)-b));
 end
end
 cdif=lo-cl;
  % check uniquness
 nc=0;
 po=zeros(1,length(s));
 for vv=2:length(cdif)
     if cdif(vv-1)>0 && cdif(vv)<=0
         nc=nc+1;
         po(nc)=vv;
     end
 end
 if nc==1
     Ccut(t) = s(po(nc));
 elseif nc==0
     Ccut(t) = 1.1;
 else
     disp('multi-cut for cp');
     t
     test=po(po>0);
     length(test)
     for qq=1:length(test)
         [lo(test(qq)) cl(test(qq))];
     end
     Ccut(t) = s(po(1));
 end
 for vv=1:length(s)
     if s(vv)<=Ccut(t)
         Vc_3(vv,t)=lo(vv);
     else
         Vc_3(vv,t)=cl(vv);
     end
 end
 % m's payoff at stage 3
 for vv=1:length(s)
     if s(vv)<=Ccut(t)
         Vm_3(vv,t)=Vm_4(vv,t);
     else
         Vm_3(vv,t)=-(b+k(t))*min(s(vv),1)+k(t);
     end 
 end
 
 % cp's decision at stage 2. 
  for jj=1:length(s)
      s2=s(jj);
     % e=h
      s3_h=s2+H;
      p3_h=min(jj+h,1000);
      s3_l=s2+L;
      p3_l=min(jj+l,1000);
      Eh = ph*Vc_3(p3_h,t)+(1-ph)*Vc_3(p3_l,t);
     % e=l
      El = Vc_3(p3_l,t);
     % decision 
      if Eh>=El
          Vc_2(jj,t)=Eh;
          Eff(jj,t)=1;
      else
          Vc_2(jj,t)=El;
          Eff(jj,t)=0;
      end
  end
  
  % m's payoff at stage 2
   for vv=1:length(s)
       
     if Eff(vv,t)==1
         Vm_2(vv,t)=ph*Vm_3(min(vv+h,1000),t)+(1-ph)*Vm_3(min(vv+l,1000),t);
     else
         Vm_2(vv,t)=Vm_3(min(vv+l,1000),t);
     end 
   end
end   
 
% Find critial points 
  % t_hat
AA = Mcut-Ccut;
BB = zeros(length(AA),1);
AA(AA<=0)=-1;
AA(AA>0)=1;
t_hat = find(AA>0,1)
  % t_m
t_m = find(Mcut>1,1)

% graph
 % conflict thresholds
CC = Mcut;
d=(Mcut(t_hat)-Mcut(1))/(t_hat-1);
V1=Mcut(1):d:Mcut(t_hat);
CC(1:t_hat)=V1;
for i=1:length(Mcut) 
    if Mcut(i)-Ccut(i)>0.045 || Mcut(i)>=1
        CC(i)=Ccut(i);
    end
end
%
  % power increase
PI = (ph*(H-L)-L)*(1:T)'+S_bar;

for i=1:T
    if PI(i)+0.01>=CC(i) && i<T
        PI(i+1)=PI(i)+L;
    end
end
for i=1:T
    if i>=t_hat && i<T
        PI(i+1)=PI(i)+ph*(H-L)-L;
    end
end
QQ = CC-PI;
End = find(QQ<=0,1);

plot(1:End,CC(1:End));
hold on
plot(1:t_hat,CC(1:t_hat),'r')
plot(1:End-1,PI(1:End-1),'g');
hold off

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Figure A.3: slow health deteroriate, conflict
clear

T=1000;
b=10;
del=0.7;
R=10;
w=0.05;
eta = 0.7;
h = 10;
H = h*0.001;
l = 1;
L = l*0.001;
r = 0.1;
p = 0:0.001:1;
s=0:0.001:1;
q=1-p;
ph=0.5;
 S_bar=0.01;

k=zeros(1,T);
for i=1:T
    AA=0;
    for j=i+1:T               
        A=del^(j-i)*prod(q(i+1:j));
        AA=AA+A;
    end
    k(i) = r+r*AA;
end


Vm_4=zeros(length(s),T);
Vm_3=zeros(length(s),T);
Vm_2=zeros(length(s),T);
Vm_1=zeros(length(s),T);
Mcut = zeros(T,1);
Vc_4=zeros(length(s),T);
Vc_3=zeros(length(s),T);
Vc_2=zeros(length(s),T);
Vc_1=zeros(length(s),T);
Ccut = zeros(T,1);
Eff=zeros(length(s),T);
for t=1000:-1:1
 % ruler's decision at stage 4
  % if keeps,
 keep=zeros(1,length(s));
 for jj=1:length(s)     
     s4 = s(jj);
     if t<1000
         keep(jj) = r+del*(p(t+1)*min(s4+w,1)*eta*R+(1-p(t+1))*Vm_2(jj,t+1));
     else
         keep(jj) = r+del*(p(t+1)*min(s4+w,1)*eta*R);
     end
 end
 strip = k(t);  
 dif=keep-strip;
 % check uniquness
 nc=0;
 po=zeros(1,length(s));
 for vv=2:length(dif)
     if dif(vv-1)>0 && dif(vv)<=0
         nc=nc+1;
         po(nc)=vv;
     end
 end
 if nc==1
     Mcut(t) = s(po(nc));
 elseif nc==0
     Mcut(t) = 1.1;
 end
 for vv=1:length(s)
     if s(vv)<Mcut(t)
         Vm_4(vv,t)=keep(vv);
     else
         Vm_4(vv,t)=strip;
     end
 end
% cp's payoff at stage 4
for jj=1:length(s)
    s4=s(jj);
     if s4>Mcut(t)
         Vc_4(jj,t)=-b;
     else
         if t<1000
             Vc_4(jj,t)=del*(p(t+1)*(min(s4+w,1)*(R+b)-b)+(1-p(t+1))*Vc_2(jj,t+1));
         else
             Vc_4(jj,t)=del*(p(t+1)*(min(s4+w,1)*(R+b)-b));
         end
     end
 end

%cp's decision at stage 3
lo = zeros(1,length(s));
cl = zeros(1,length(s));
for jj=1:length(s)
 % challenge
 s3=s(jj);
 cl(jj) = min(s3,1)*(R+b)-b;
 % loyal
 if t<1000
     lo(jj) = Vc_4(jj,t);
 else
     lo(jj) = del*(p(t+1)*(min(s3+w,1)*(R+b)-b));
 end
end
 cdif=lo-cl;
  % check uniquness
 nc=0;
 po=zeros(1,length(s));
 for vv=2:length(cdif)
     if cdif(vv-1)>0 && cdif(vv)<=0
         nc=nc+1;
         po(nc)=vv;
     end
 end
 if nc==1
     Ccut(t) = s(po(nc));
 elseif nc==0
     Ccut(t) = 1.1;
 else
     test=po(po>0);
     for qq=1:length(test)
         [lo(test(qq)) cl(test(qq))];
     end
     Ccut(t) = s(po(1));
 end
 for vv=1:length(s)
     if s(vv)<=Ccut(t)
         Vc_3(vv,t)=lo(vv);
     else
         Vc_3(vv,t)=cl(vv);
     end
 end
 % m's payoff at stage 3
 for vv=1:length(s)
     if s(vv)<=Ccut(t)
         Vm_3(vv,t)=Vm_4(vv,t);
     else
         Vm_3(vv,t)=-(b+k(t))*min(s(vv),1)+k(t);
     end 
 end
 
 % cp's decision at stage 2. 
  for jj=1:length(s)
      s2=s(jj);
     % e=h
      s3_h=s2+H;
      p3_h=min(jj+h,1000);
      s3_l=s2+L;
      p3_l=min(jj+l,1000);
      Eh = ph*Vc_3(p3_h,t)+(1-ph)*Vc_3(p3_l,t);
     % e=l
      El = Vc_3(p3_l,t);
     % decision 
      if Eh>=El
          Vc_2(jj,t)=Eh;
          Eff(jj,t)=1;
      else
          Vc_2(jj,t)=El;
          Eff(jj,t)=0;
      end
  end
  
  % m's payoff at stage 2
   for vv=1:length(s)
       
     if Eff(vv,t)==1
         Vm_2(vv,t)=ph*Vm_3(min(vv+h,1000),t)+(1-ph)*Vm_3(min(vv+l,1000),t);
     else
         Vm_2(vv,t)=Vm_3(min(vv+l,1000),t);
     end 
   end
end  

% Find critial points 
  % t_hat
AA = Mcut-Ccut;
BB = zeros(length(AA),1);
AA(AA<=0)=-1;
AA(AA>0)=1;
t_hat = find(AA>0,1)

  % t_m
t_m = find(Mcut>1,1)  

% graph
 % conflict thresholds

CC = Mcut;
d=(Mcut(t_hat)-Mcut(1))/(t_hat-1);
V1=Mcut(1):d:Mcut(t_hat);
CC(1:t_hat)=V1;
for i=1:length(Mcut) 
    if Mcut(i)-Ccut(i)>0.045 || Mcut(i)>=1
        CC(i)=Ccut(i);
    end
end

  % power increase
PI = (ph*(H-L)-L)*(1:T)'+S_bar;
for i=1:T
    if PI(i)+0.01>=CC(i) && i<T
        PI(i+1)=PI(i)+L; 
    end
end
for i=1:T
    if PI(i)+ph*(H-L)-L>=CC(i) && i<T
        PI(i+1)=PI(i)+ph*(H-L)-L;
    end
end
QQ = CC-PI;
End = find(QQ<=0,1);

plot(1:End-1,CC(1:End-1));
hold on
plot(1:t_hat,CC(1:t_hat),'r')
plot(1:End-1,PI(1:End-1),'g');
hold off

diary off